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Starting from a Cluster Variational Method, and inspired by the correctness of the param- 
agnetic Ansatz (at high temperatures in general, and at any temperature in the 2D Edwards- 
Anderson model) we propose a novel message passing algorithm — the Dual algorithm — 
to estimate the marginal probabilities of spin glasses on finite dimensional lattices. We show 
that in a wide range of temperatures our algorithm compares very well with Monte Carlo 
simulations, with the Double Loop algorithm and with exact calculation of the ground state 
of 2D systems with bimodal and Gaussian interactions. Moreover it is usually 100 times 
faster than other provably convergent methods, as the Double Loop algorithm. 
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I. 



INTRODUCTION 



Inference problems are common to almost any scientific discipline. Often an inference problem 
can be recast in that of computing some marginal probability on a small subset of variables, given 
a joint probability distribution over a large number N of variables, {si, . . . ,S]\f} = a. A typical 
example is the computation of the expectation value of a variable (sj) = Yls ^iPii^i) where Pi{si) 
is the single-variable marginal probability defined as 



Clearly the computation of the sum in the r.h.s. is as difficult as the computation of the partition 
function Z = which is the main subject of statistical mechanics. In the general (and the 

most interesting) case, this problem can not be solved exactly in a time growing sub-exponentially 
with the size A'^. We are deemed to use some kind of approximation in order to compute marginals 
in a time growing linearly in N . The approximation schemes used so far are mainly adopted from 
the field of statistical mechanics [H, where mean-field-like approximations are standard and well 
controlled tools for approximating the free-energy, —f5f = In Z. 

For non-disordered models, like the Ising ferromagnetic model, these approximations work quite 
well and provide good results 0]. However, much less is known for models with disorder, and for 
this reason we will focus on spin glasses in the present paper. In disordered models one usually deals 
with an ensemble of problems (e.g. in spin glasses each sample has its own couplings, randomly 
chosen from a given distribution), and the results obtained by the statistical mechanics tools refer 
to average quantities, i.e. those of the typical samples. In other words one is not concerned with 
the behavior of a specific sample, but rather one looks at the whole ensemble. On the contrary 
when doing inference one is interested on the properties of a single specific problem and thus the 
above approximation schemes (based on statistical mechanics mean-field approaches) need to be 
converted into an algorithm that can be run on such a specific sample. Computing marginals on a 
given sample clearly gives more information than computing averages over the ensemble. 

To our knowledge, effective (i.e. linear time in N) algorithms for computing marginals can be 
essentially divided in two broad classes: stochastic local search algorithms, that roughly sample 
the configurational space according to -P(c), and algorithms based on some kind of mean-field 
approximation. The former are exact on the long run, but the latter can be much more useful if an 
approximated answer is required in a short time. Unfortunately the latter also have some additional 
drawback due to the mean-field nature of the underlying approximation, e.g. the appearance of 
spurious phase transitions, that may prevent the proper convergence of the algorithm. 

One more reason why this latter class of algorithms has a limited scope of application is that the 
convergence of the algorithm may strongly depend on the presence of short loops in the network of 
variables interactions. In this sense the successful application of one of these algorithms to models 
defined on regular lattices (which have many short loops) would be a major achievement. 

In this paper we introduce a fast algorithm for computing marginals in 2D and 3D spin glass 
models 0] defined by the Hamiltonian (further details are given below) 



The first non trivial mean-field approximation for the above model corresponds to the Bethe- 
Peierls approximation scheme and the Belief Propagation (BP) algorithm Q. Unfortunately when 
BP is run on a given spin glass sample defined on a -D dimensional lattice, it seems to provide 
exactly the same output as if it were run on a random regular graph with fixed degree 2D: that 
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is, for T > TBethe it converges to a solution with all null local marginals ((sj) = 0), while for 
T < Teethe it does not converge to a fixed point. 

The next step in the series of mean-field approximations (also known as Kikuchi approxima- 
tions or 'cluster variation method') is to consider joint probability distributions of the four spins 
belonging to the same plaquette [1, S]- Under this approximation an algorithm has been derived 
which is called Generalized Belief Propagation (GBP) To our knowledge, this algorithm has 
been applied to 2D spin glasses, but only in presence of an external magnetic field, which is known 
to improve the convergence properties of GBP. Our experience says that, running plain GBP on 
a generic sample of a 2D spin glass without external field, a fixed point is reached only for high 
enough temperatures, well above the interesting region. The lack of convergence of GBP (and other 
similar message-passing algorithms, MPA) is a well known problem, whose solution is in general far 
from being understood. For this reason, a new class of algorithms has been recently introduced [3], 
which provably converge to a fixed point: these algorithms use a double loop iterative procedure 
(to be compared to the single loop in GBP) and for this reason they are usually quite slow. 

The algorithm we are going to introduce is as fast as BP (which is the fastest algorithm pre- 
senty available), and converges in a wider range of temperatures than BP. Moreover the marginal 
probabilities provided by our algorithm are as accurate as those which can be obtained by Double 
Loop algorithms in a much larger running time. 

Our algorithm works in the absence of external field, that is in the situation where MPA have 
more problems in converging to a fixed point. In this sense it is a very important extension to 
presently available inference algorithms. 

The rest of the work is organized as follows. First, in section HI] we derive the GBP equations 
for the 2D Edward Anderson model. Then, in section III Al we rewrite these equations in terms 
of fields, a notation that has a nicer physical interpretation and that we are going to use in the 
rest of the work. Section IIIII presents the novel algorithm, inspired by the paramagnetic Ansatz 
to the GBP equations. In section IIVI we show the results of running this algorithm on the 2D 
Edward-Anderson model. There we compare its performance with Monte Carlo simulations, with 
the Double Loop algorithm and with exact calculation of the ground state of systems with bimodal 
and Gaussian interactions. Then, in section |V] we generalized our message passing equation to 
general dimensions and present some results for the 3D Edward-Anderson model. Finally, some 
conclusions are drawn in section IVIi 



II. GENERALIZED BELIEF PROPAGATION ON THE 2D EA MODEL 

Here we present the GBP equations for the Edwards-Anderson (EA) model on a 2D square 
lattice, and we refer the reader to 0, 0] for a more general introduction. In our case (as well as in 
many other cases) GBP is equivalent to Kikuchi's approximation, known as the Cluster Variational 
Method (CVM) [5|. We will try a presentation as physical as possible. 

Consider the 2D EA model consisting of a set a = {si, . . . , sj\f} N Ising spins Sj = ±1 located 
at the nodes of a 2D square lattice, interacting with a Hamiltonian 

^(o-) = - ^ JijSiSj , (2) 

{id) 

where the sum runs over all couples of neighboring spins (first neighbors on the lattice). The Jij 
are the coupling constants between spins and are supposed to be fixed for any given instance of 
the model. If the interactions are not random variables, i.e. Jij = J, then the 2D ferromagnet is 
recovered. We will focus on the two most common disorder distributions: bimodal interactions, 
P(J) = l5{J - 1) + ^6{J + 1), and Gaussian interactions P{J) = exp(-j2/2)/\/27r. 
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The statistical mechanics of the EA model, at a given temperature T = 1//3, is given by the 
Gibbs-Boltzmann distribution 

P(^) = — . 

The direct computation of the partition function 

(7 

or any marginal distribution p{si,Sj) = sj ^i^) is a time consuming task, unattainable in 

practice, since it involves the addition of 2^ terms, and therefore an approximation is requirecl^. 

The idea of the Region Graph Approximation to Free Energy [6,] is to replace the real dis- 
tribution P{cr) by a reduced set of its marginals. The hierarchy of approximations is given by 
the size of such marginals, starting with the set of all single spins marginals Pi{si) [mean- field] , 
then following to all neighboring sites marginals Pij{si,Sj) [Bethe], then to all square plaquettes 
marginals pijki{si, Sj, Sk, si), and so on. Since the only way of knowing such marginals exactly is 
the unattainable computation of Z, the method pretends to approximate them by a set of beliefs 
bi{si), bij{si, Sj), etc. obtained from the minimization of a region based free energy. In the Region 
Graph Approximation to Free Energy, a set of regions, i.e. sets of variables and their interactions, is 
defined, and a free energy is written in terms of the beliefs at each region. The Cluster Variational 
Method does a similar job, but instead of starting from an arbitrary choice of regions, it starts by 
defining the set of largest regions, and smaller regions are defined recursively by the intersections 
of bigger regions. In this sense, CVM is a specific choice of region graph approximation to the free 
energy. 

For the 2D EA model, we will consider the expansion of the free energy in terms of the marginals 
at three levels of regions: single sites (or spins), links, and plaquettes. By plaquettes we mean the 
square basic cell of the 2D lattice. This choice of regions corresponds to the CVM having the 
square plaquettes as biggest regions. The free energy of the system is therefore written as 



Y^crFr 



R 

where R runs over all regions considered, and the free energy in a particular region depends on the 
marginals at that level 6_R(cri?): 

PFr = bR{aR)(3ER{aR) + hR^aR) log 6ij(crR) . 

The symbol aR refers to the set of spins in region R, while Er is the energy contribution in that re- 
gion. The counting numbers cr (also Moebius coefficients) are needed to ensure that bigger regions 
do not over count the contribution in free energy of smaller regions, and follow the prescription 

CR = l- Y Cri , (3) 

where R' is any region containing completely region i?, as e.g. a plaquette containing a link or a 
link containing a site. In the case of the 2D lattice, the biggest regions are the square plaquettes, 



^ The 2D case is special: indeed, thanks to the small genus topology, the partition function Z can be computed 
efficiently. Ifowever we are interested in developing an algorithm for the general case, and we will not make use of 
this peculiarity. 
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and therefore Cpiaq = 1, while the hnks regions have cnnk = 1 — 2cpiaq = — 1 (as each of them is 
contained in 2 plaquettes regions) , and finally the spins regions have Cgite = 1 — 4cpiaq — 4ciink = 1 
(as each spin belongs to 4 links and 4 plaquettes). So the actual approximation for the EA model 
on a 2D square lattice is 



/3F = Y,J2^V M log ^^^l^] Plaquett. 



es 



EE^^("^)l°g r^^^'l ,1 Links (4) 
lT^ exp[-/?SL(cJL)] 



biisi 



+ E E ^^(^^) r\i^f V Sites 

W exp[-^£;i(si)J 



where the sums run over all plaquettes, links and sites respectively. Please note that we are using 
the following notation for region indices: lower-case for sites, upper-case for links and upper-case 
calligraphic for plaquettes. The energy term Ei{si) in the sites contribution is only relevant when 
an external field acts over spins, and can be taken as zero in our case, since no external field is 
considered. Notice that whenever the interactions are included in more than one region (in our 
case are included in Link and Plaquette regions), the counting numbers guarantee that the exact 
thermodynamical energy U = '^(j P{(y)T~i{a) is obtained when the beliefs are the exact marginals 
of the Boltzmann distribution. On the other hand, the entropy contribution is intrinsically ap- 
proximated, since the cutoff in the regions sizes imposes a certain kind of factorization of P{cf) in 
terms of its marginals (see for an explanation of the region graph approximation in terms of 
cumulants expansions of the entropy) . 

The next step in the method, is to compute the beliefs from the minimization condition of 
the free energy. However, an unrestricted minimization will generally produce inconsistent solu- 
tions, since the beliefs (marginals) are not independent, as they are related by the marginalization 
conditions 



Hsi) = Eai\, blioL) = E., bL{Si,Sj) , 

blicTL) = bLiSi,Sj) = T.aj,\r. ^'pM = Ssfc.s, ^risi, 



(5) 



where ctl = {sj, sj} and a-p = {si, Sj, s^, s/}. In order to minimize under the constraints in Eq. ^ 
and under the normalization condition for each belief, a set of Lagrange multipliers should be 
added to the free energy in Eq. There are different ways of choosing the Lagrange multipliers 
0], and each of them will produce a different set of self consistency equations. We choose the so 
called Parent to Child scheme (see section XX in 0]), in which constraints in Eq. ([5]) are imposed 
by two sets of Lagrange multipliers: /XL_>.j(sj) relating the belief at link L to that at site i, and 
i''p^L{aL) relating the one at plaquette V to the one at link L. 

With constraints ([5]) enforced by Lagrange multipliers, the free energy stationary conditions for 
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the beliefs are the fohowing: 



biisi 



^exp ( -/3Si(si) - y'/XL^i(si) ) , 



Zl 



exp 



exp 



V 



2 3 



VDL 



iCL L'Di 



4 2 



LcrV'DL 
V'+V 



L(tV 



(6) 



J 



where the notation L D i refers to all links containing site i and V D L to all plaquettes containing 
link L. The upper indices in the sums are written just to help understanding how many terms 
are in each sum for the 2D case. The precise meaning of the indices in each summation can 
be understood from the graphical representation in Figure [TJ Lagrange multipliers are shown as 
arrows going from parent regions to children regions: simple arrows correspond to and triple 

arrows to f-p-^L- Let us consider, for instance, the belief in a link region biicTL), depicted in the 
central picture of Fig. [TJ the sum of the two Lagrange multipliers u-p^L^ai,) corresponds to the 
triple arrows from plaquettes on the left and right of the link L, while the double sum over the 
three fiL^i^Si) and the three HL^j{sj) correspond to the three arrows acting over the two spins. 



Zl 




FIG. 1. Schematic representation of belief equations (jl]). Lagrange multipliers are depicted as arrows, going 
from parent regions to children regions. 



In Eq.(l6]), the Z/j are normalization constants, and the terms E-p^a-p) = Ep{si, sj, Sk, si) = 
-{JijSiSj + JjkSjSk + JkiSkSk + JiiSiSi) and Ei^aL) = EL{si, sj) = - JijSiSj are the corresponding 
energies in plaquettes and links, and are represented in Fig. [Uby bold lines (interactions) between 
circles (spins). In our case Ei{si) is zero since no external field is acting on the spins. 

The Lagrange multipliers are fixed by the constraints they were supposed to enforce, Eq. ([5]), 
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and they must satisfy the following set of self-consistency equations: 

2 3 



exp 



exp[^ - iyp^L{si, Sj) - fiD^iisi) - fiu^j{sj)j 



L'Dj 



(7) 



E 



exp 



L'eVV'DL' 
L'j^L V'+V 



L'Dl 



Again, to help understanding these equations, we provide in Fig. [2] their graphical representation. 
Note that there is one of these equations for every pair of Link-Site and every pair of Plaquette-Link 
in the graph. With Ep\j^ we refer to interactions in plaquette P that are not in link L. 



© 







•D 



FIG. 2. Message passing equations ([7]), shown schematically. Messages are depicted as arrows, going from 
parent regions to children regions. On any link Jij, represented as bold lines between spins (circles), a 
Boltzmann factor e^'^'J**^^ exists. Dark circles represent spins to be traced over. Messages from plaquettes 
to links i^p^L^Si, Sj) are represented by a triple arrow, because they can be written in terms of three 
parameters U, Ui and Uj, defining the correlation {siSj) and magnetizations (si) and (sj), respectively. 



For each link L in the 2D lattice, there are two link-to-site multipliers, m^i{si) and fiL^j{sj). 
For each plaquette there are four plaquette-to-link multipliers vp^L^Si, Sj), corresponding to the 
four links contained inside the plaquette. Let be the number of spins in the lattice, there are 2N 
links and plaquettes. So the originally intractable problem of computing marginals, has been 
replaced by the problem of solving a set of AN + AN coupled equations for the Lagrange multipliers 
as those in Eq. ([7|). Once these equations are solved, the approximation for the marginals is 
obtained from Eq. ([6]) for the beliefs, and all thermodynamic quantities are derived from them as 
in Eq. (jH). 

Minimizing a region graph approximation to free energy, as that in Eq. with constraints 
Eq. ([5]), or equivalently solving the set of self-consistent equations in Eq. d?]), is still a non trivial 
task. Let us consider two ways of doing it. The first method is the "direct" minimization of the 
constrained free energy, using a Double Loop algorithm 0]. This method is quite solid, since it 
guarantees convergence to an extremal point of the constrained free energy, but it may be very slow 
to converge. The second method, which is generally faster but is not guaranteed to converge, is the 
family of the so called Message-Passing algorithm (MPA), in which the Lagrange multipliers are 
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interpreted as messages v-p^i^ai) going from plaquettes to links, and messages ^L-fi{si) from links 
to sites. Self consistency equations ([7]) can be viewed as the update rules for the messages in the 
left hand side, in terms of those in the right hand side. A random order updating of the messages 
in the graph by Eq. ([7]) (message passing) can reach a fixed point solution, and therefore, to an 
extremal point of the constrained free energy 0]. Next, we show explicitly how the message-passing 
equations looks like in terms of fields. 



A. Prom multipliers to fields 

A particularly useful way of representing the multipliers (messages), with a nice physical inter- 
pretation, is the one used in which we adopt here. In full generality [3, 0], these multipliers 
can be written in terms of effective fields 

IJ.L^i{Si) = /3 UL^i Si (8) 
J^V^LiSi, Sj) = (3 (U-p^L Si Sj + U-p^i Si + U-p^j Sj) (9) 



In particular, the field u corresponds to the cavity field in the Bethe approximation [6[. Using 
Lagrange multipliers, messages or fields, is essentially equivalent. We will often refer to fields 
as n-messages to emphasize their role in a message-passing algorithm, and we will refer to self- 
consistency equations d?]) as the message-passing equations. 

This parametrization of the multipliers has proved useful to other endeavors, like the extension 
of the replica theory to general region graph approximations Here, all the relevant information 
in the Lagrange multipliers is translated to "effective fields" u and {U,Ua,ui,). Notice that in 
this representation every single field u corresponds to an arrow in the schematic messages-passing 
equations in Figure [21 In particular, the messages going from plaquettes to links are characterized 
by three fields {U,Ua,Ub), and the field U acts as an effective interaction term, that adds directly 
to the energy terms in the Boltzmann factor. For instance, the first message-passing Eq. ([7]) is 



exp 



/SuL^iS^^ = ^ exp 



l3[{u'p^i + Uc^i) Si + [U-p^L + Uc^L + Jij) SiSj + 
{u-p^j + Uc^j + UA^j + UB^j + Uu^j) Sj 



(10) 



where the indices refer to the notation used in Fig. [2] and Jij is the interaction coupling constant 
between spins Si and Sj. This equation naturally defines the updating rule for the message UL^f. 

UL^i = u{up^i + Uc-^i, Up^L + Uc^L + Jij, Up^j + Uc-^j + UA^j + UB^j + Uu^j) , (H) 

where 

.(n,^,/.)=n+-log-^^^^^ 



Note that the usual cavity equation for fields in the Bethe approximation 10(] is recovered if all 
contributions from plaquettes V and £ are set to zero. 

Working in a similar way for the second equation in ([7]) we end up with the updating rule for 
the message {Upi^L,up^i,up^j) sent from any given plaquette region V to one of its children 
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links L (see right picture in Fig. [2]) 



Uv^L 




(12) 



where 



K{si,Sj) = ^ exp (iUUu^u + Jjk)sjSk + (U-ji-^r + Jki)skSi + (Uv-^d + Jii)siSi + 



{uu^k + UC^k + UE~^k + UR-^k)sk + {UR^I + UF^l + UG^i + UD^l)si 



Equations (jlip and ()12p are equivalent to equations in ([7]), once multipliers (messages) are 
parametrized in terms of fields. For instance, note that the fi multipliers in the left hand side of 
second equation in ^ appear now subtracted in the right hand side of Eq. ([12]) . 

The field notation is more comprehensible and has a clear physical meaning. Each plaquette 
V is telling its children links L that they should add an effective interaction term U-p^i to the 
real interaction Jjj, due to the fact that spins Si and sj are also interacting through the other 
three links in the plaquette V. Fields u act like magnetic fields upon spins, and the complete 
i^v^Lisi, Sj)— message is characterized by the triplet {U-p^L,u-p^i,u-p^j), and will be referred from 
now on as ?7nn— message. Furthermore, it is clear that some fields enter directly in the message- 
passing equations like up^i and uc^i in Eq. ([TT]) and UT>^i and uu^j in Eq. (fT2]l . Also note that 
since our model has no external field, the fields u break the symmetry of the original Hamiltonian 
whenever they are non zero. For instance, in the ferromagnetic model, when all Jij = J, these 
fields are zero at high temperature and become non zero at Kikuchi's critical temperature T = 1.42 
0], implying a spontaneous magnetization in the ferromagnet. 

III. THE DUAL APPROXIMATION FOR THE PARAMAGNETIC PHASE 

Unfortunately, the iterative message-passing algorithm for solving the GBP equations (llip and 
(jl2p often does not converge on finite dimensional lattices. While this is expected if long range 
correlations are present, it is rather disappointing that it happens also in the paramagnetic phase, 
where one would like to find easily the solution to the model. Here we are going to focus only on the 
paramagnetic phase, and propose an improved solving algorithm based on physical assumptions. 

In the paramagnetic phase of any spin model defined by the Hamiltonian in Eq. ([2]), that is 
with no external field, variables have no polarization or magnetization: this in turn implies that in 
the solution all «— fields must be zero, and only [/—fields should be fixed self-consistently to non 
zero values. 

This paramagnetic solution has some interesting properties. First, it is always a solution of 
the GBP equations, since Eq. (jlip and Eq. (|12p are self-consistent with all u = 0. This means 
that starting from unbiased messages (all u = 0) the iterative GBP algorithm keeps this property. 
Second, the paramagnetic Ansatz is correct, from the GBP perspective, at least at high enough 
temperatures, meaning that even if we start with biased messages {u ^ 0), the iterative algorithm 
converges to all u = at high temperatures. And last, but not least, the well studied physical 
behavior for the 2D EA model with zero- mean random interactions Jij, is expected to remain 
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always paramagnetic, i.e. to have no transition to a spin glass phase at any finite temperature ll| 
Therefore, the Ansatz u = is both physically plausible and algorithmically desirable. 



® 
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FIG. 3. Message passing of correlation messages in the dual approximation. In the right hand side the trace 
is taken over the black spins. 



Under the paramagnetic Ansatz, which we shall also call Dual approximation for a reason to be 
explained soon, the message-passing equation (llip is irrelevant, as it is always satisfied given that 
u(0, [/, 0) = 0, while Eq. now turns into (see Fig. [3]) 

1 



— arctanh 



tanh l3{Uu^u + J jk)tanhP{Un^R + JM)ta.nhl3{Uv^D + Jh) ■ (13) 



The only relevant messages now are those associated to the multipliers ^^^^^{si, Sj) = jiU-p^i Si Sj, 
and they will be refereed to as [/—messages. Eq. (fT3]) can be interpreted as a correlation message- 
passing equation, giving the new interaction field U that a certain link shall experience as a 
consequence of the correlations transmitted around the plaquette. The belief Eq. ([6]) also simplify. 
Obviously 6(sj) = 0.5 for every spin in the graph, and the link and plaquette beliefs are 

biisi, Sj) = ^e^iUc^L+Uv^i.+J,,)s,si ^ ^j4) 

b'p{Si,S-,Sk,Si) = -}—QPi^C^L+Jij)siSj+l3{Uu^u+Jjk)sjSk+P(UTz^R+Jki)skSi+li{U'D^D+Jii)siSt 

Z'p 

The Dual algorithm we are proposing to study the paramagnetic phase of the EA model, is a 
standard message passing algorithm for the [/-messages, which works as follows. 
1: Start with all [/-messages null 
2; repeat 

3; Choose randomly one plaquette V and one of its children links L 

4: Update the field U-p^L according to Eq. (fT3]) as in Fig. [3] 

5: until The last change for any [/-message is less than e (we use typically e = 10^^") 

6: return The beliefs bi^Si, sj) defined in Eq. (jl4p for every pair of neighboring spins 

Some damping factor 7 G [0, 1) can be added in the update step U-p^i = jU-p-^L + (1 — 7)[/ in 
order to help convergence. 
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A. Mapping to the dual model 



It is worth noticing that Eq. ()13p is nothing but the BP equation for the corresponding dual 
model (hence the name of the algorithm). 

The dual model has a binary variable Xij = SiSj on every link of the original model, and the 
original coupling constants play now the role of an external polarizing (eventually random) field 



^dual(^) = -"^Jj 



ijXij . 



This Hamiltonian looks like the sum of independent variables, but this is not the case. The dual 
variables Xij = ±1 must satisfy a constraint for each cycle (or closed path) in the original graph, 
enforcing that their product along the cycle must be equal to 1. On a regular lattice any closed 
path can be expressed in terms of elementary cycles of 4 links (the plaquettes) and so it is enough 
to enforce the constraint on every plaquette: XijXjkXtixu = 1. The Gibbs-Boltzmann probability 
distribution for the dual model is then given by 

P^g) = -Q (15) 

where the product runs over all elementary plaquette. 

The model described by the probability measure in Eq. ()15p can be viewed as a constraint 
satisfaction problem with a non uniform prior (given by e~^^duai(^))_ \i straightforward to derive 
the BP equations for such a problem. Indeed by defining the marginal for the variable Xij on link L 
in the presence of the solely neighboring plaquette V as {l + Xij tanh fUJ-p^i) /2 oc exp{f5U-p^LXij) , 
the BP equations read 

^ (l + Xij tanh /3C/-p^l) oc 

■^jk ^^kli-^li 

^ (l + Xjk tanh (3{Uu-^u + Jjk)) {I + xh tanh /3{Un^B^ + Jki)) (l + xu tanh P{Uv~^d + Ju)) = 

^jki-^kl i^li ■ 
•^jk'^kl-^li — "^ij 

1 + Xij tanh f3{Uu-yu + Jjk) tanh f3{UTi^R + Jki) tanh f3{Uv^D + Jh) ■ (16) 

In the second summation the terms containing one or two x variables sums to zero, while the other 
two terms are those written in the last expression. Equating the first and the last expressions, this 
equation is manifestly equal to Eq. (fT3|) . 



B. Average Case Solution 

GBP in general, and the Dual approximation in particular, are methods for the study of the 
thermodynamical properties of a given problem. However, in the limit of large systems {N — )• oo, 
thermodynamical limit), we expect a typical behavior to arise. This is the so called self-averaging 
property of disordered systems. By typical we mean that almost every realization of the interactions 
Jij will result in a system whose thermodynamical properties (free energy, energy, entropy) are very 
close to the average value. 

Normally, in disordered systems, we cope with the N ^ oo limit and with the average over 
the random Jij by the replica method. The application of the replica trick to regions graph 
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approximations is a challenging task [9|. However, we can still grasp the average case behavior 
with a cavity average case solution of the dual message-passing equations, at the price of neglecting 
the local structure of the graph (beyond plaquettes). 

The idea is to represent the set of ?7-messages flowing in any given graph, by a population of 
messages Q{U). Then the message-passing Eq. ([13]) is used to obtain such population in a self- 
consistent way. More precisely, in every iteration three messages Ui,U2,Us are randomly drawn 
from the population Q{U) and a new message Uq = U{Ui, U2, U3) is computed by Eq. (fT3l) using 
three couplings randomly selected from P{J). The obtained message Uq is put back into the 
population, and the iteration is repeated many times, until the population stabilizes. 

Once we have the self consistent population of messages, we can compute the average energy 

^Ave = tanh^ ( + Ui + U2))Q(u,),Q{U2),P{J^j) (1^) 

by a random sampling of the population and of the interactions. The average case solution is 
supposed to be very good whenever the network of interactions has no or few short loops. This is 
not the case in any finite dimensional lattice, since there the short loops (plaquettes) are abundant. 
Nonetheless, the average case solution gives a reasonably good approximation to the single instance 
results in 2D and 3D, as shown in the next Section. 

IV. RESULTS ON 2D EA MODEL 

Message-passing algorithms work fine in the high temperature regime (T > Tc) of models de- 
fined on random topologies: this is the reason why these methods have been successfully applied 
in random constraint satisfaction problems, like random-SAT or random-Coloring How- 
ever, when used on regular finite-dimensional lattices, they can experience difficulties even in the 
paramagnetic phase, because the presence of short loops spoils message-passing convergence. 

It is well known that on a random graph of fixed degree (connectivity) c = 4 the cavity approx- 
imation gives a paramagnetic result above Teethe — 1-52 (i.e. f^sethe — 0.66) with all cavity fields 
Ui = 0. Below the Bethe critical temperature, this solution becomes unstable to perturbations, 
and we expect many solutions to appear with non trivial messages Ui 7^ 0. The presence of many 
solutions in the messages passing equations is connected to the existence of many thermodynami- 
cal states in the Gibbs-Boltzmann measure, or, equivalently, to the presence of replica symmetry 
breaking. The appearance of such a spin glass phase is also responsible for the lack of convergence 
of message-passing equations, since the intrinsic locality of the message-passing equations fails to 
coordinate distant regions of the graph (which are now long-range correlated). As a consequence, 
the application of BP to the 2D EA model (that also has fixed degree c = 4) still finds the paramag- 
netic phase at high temperatures, but below Teethe , the Bethe instability takes the message-passing 
iteration away from the u = solution and does not allow the messages to convergence to a fixed 
point {i.e. the algorithm wanders forever). In Figure [4] we show the convergence probability for 
the BP message-passing equations in the 2D EA model. 

On the other hand, a straightforward GBP Parent-to-Child implementation does not fully over- 
come this problem. At high temperatures, the Parent-to-Child equations converge to a param- 
agnetic solution with all u = and non trivial U ^ 0, which turns out to be the same solution 
found by our Dual algorithm. When going down in temperature, the convergence properties of the 
algorithm worsen, and are sensitive to tricks like damping and bounding in the fields. A thorough 
discussion of these properties is left for a future work, but let us summarize that typically the 
algorithm stop converging at low temperatures, somewhere below Teethe; as shown in Figured 

So, in general, BP and GBP equations are not simple to use in finite dimensional systems 
at low enough temperatures: this warning was already reported in Refs. 0, 0, H]- Indeed a 
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FIG. 4. Convergence probability of BP (Bethe approximation) and GBP on a 2D square lattice, as a function 
of inverse temperature. Data points are averages over 100 systems with random bimodal interactions. System 
sizes are N = L'^ with L — 32, 128 and a damping factor 7 = 0.5 has been used in the iteration of the message- 
passing equations. The Bethe spin glass transition is expected to occur at /Seethe — 0.66 (Teethe — 1.52) for 
a random graph with the same connectivity as the 2D square lattice. Notably, that temperature also marks 
the convergence threshold for BP equations in the 2D square lattice. GBP, on the contrary, reaches lower 
temperatures, but eventually stop converging. 



different method for extremizing the constrained free energy named Double Loop algorithm 0, 
151 ] was developed to overcome such difficulties. As mentioned earlier. Double Loop guarantees 



convergence of the beliefs, on any topology, with or without short loops. Given the convergence 
problems in GBP, researchers typically resort to Double Loop algorithms to extremize region graph 
approximations to the free energy, below the Bethe critical temperature. 

In order to make a fair comparison with our Dual algorithm, we have used an optimized code 
for GBP and Double Loop algorithms: the open source LibDai library written in C-|— |- [16i] . 

The first interesting result of our work is that our Dual algorithm converges at all temperatures, 
just as Double Loop does. The reason why it converges is that there are no u-messages, so the 
Bethe instability will not affect our message-passing iteration. 

The second relevant result of our Dual algorithm, is the fact that it finds the same solution 
found by the Double Loop algorithm at all temperatures. In other words, the direct extremization 
of the region graph approximation to free energy Eq. ^ via a Double Loop algorithm finds a 
paramagnetic solution characterized by the beliefs bi{si) = 0.5 and biisi, sj) = le~^'^»j''»^j ; and the 
effective interactions Jij found by the Double Loop algorithm are exactly equal to those found with 
our Dual algorithm, Jij = Jij + U-p^i + Uc^l. This means that beliefs and correlations found by 
the two algorithms are identical: (sjSj)Doubie Loop = (■5jSj)Duai- 

The third result is that the running times of our Dual algorithm are nearly four orders of 
magnitude smaller than those required by the Double Loop implementation in libdai, at least in 
a wide range of temperatures (see figure [5]). More precisely, the convergence time of the Dual 
algorithm growth exponentially with j3 = l/T, but still, in the relevant range of temperatures 
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where the region graph approximation is a good approximation (not too low temperatures) 
running time is always roughly a factor 10^ smaller than Double Loop. 
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FIG. 5. Running times of the Double Loop algorithm 0, [IB] (libdai) and the Dual algorithm averaged over 
10 realizations of a 2D 8x8 EA model with Gaussian interactions. Generally the Double Loop algorithm 
requires a time 4 orders of magnitude larger than that used by the Dual algorithm. Three different precision 
goals where used for the Dual algorithm 10^^, 10^^", 10^^^, while the precision of the Double Loop algorithm 
is 10~^. The inset shows the behavior of the running times for both algorithm versus the system size L = \/N. 



A. Dual approximation vs Monte Carlo simulations 



The fact that our Dual algorithm provides the same results (and much faster) than the Double 
Loop algorithm is a very good news. Essentially is telling us that we are not loosing anything by 
restricting the space of possible messages, as far as the region graph approximation is concerned. 
However, the ultimate comparison for the approximation has to be done with the exact marginals 
and correlations. In figure [6] we show a comparison between the exact correlations (sjSj)pT of 
neighboring spins obtained with a Parallel Tempering (PT) Monte Carlo simulation, and the Dual 
approximation estimate for the same two-spins correlations. The coincidence between (siSj)pT and 
(■SiSj)Duai is essentially perfect at high temperatures, and it becomes weaker as the temperature 
is decreased. The reason for the discrepancies is obviously the fact that we are using an approx- 
imation in which collective behaviors of spins is accounted exactly only until the plaquette level; 
more distant correlations are approximated and these correlations become more important at low 
temperatures. 

Given such a good correspondence between the correlations under the Dual approximation and 
the true correlations, we expect a very good estimate for the energy too. In Figure [7] we show 
with points the energy under the Dual approximation and with full lines the Monte Carlo exact 
energy: the data are indeed very close. The dashed lines show the average case energy for the Dual 
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FIG. 6. Comparison between the correlations (siSj)Duai obtained by the Dual algorithm and the nearly 
exact correlations obtained by a Parallel Tempering simulation. We used a 64 x 64 EA model with Gaussian 
interactions. At each temperature the data correlation coefficient p is reported. 



approximation, Eq. (jl7p . In spite of the fact that the average case does not take into account the 
local structure of the lattice, the average case energy is quite close to the single instance one. 



B. Ground State Configuration in 2D 



The good agreement between the correlations found by the Dual algorithm, and those found in 
a Monte Carlo simulation, for the 2D EA model, compels us to push this correspondence down to 
T = 0. More precisely, using the correlations obtained by our Dual algorithm at low temperatures, 
we try to compute a ground state configuration by the following procedure. The idea is to freeze 
iteratively the relative position SjSj of those interacting spins that are more strongly correlated, 
which is done by setting Jjj — t- ±00, and re-running the Dual algorithm until convergence every 
time one pair of spins is frozen. Note that freezing the relative position of spins is equivalent to 
freeze the dual variable Xij = SiSj. The freezing procedure is very simple, but for the fact one has 
to check that frozen links must be consistent with a spin configuration. More precisely, frozen Xij 
variables must satisfy the requirement that on any closed loops the product is one, 

closed loop 

n ^^j- = 1 • (18) 

For very short loops the satisfaction of this condition is automatically induced by the Dual al- 
gorithm: for example if three links on a plaquette freeze, the fourth link is immediately frozen 
to a value satisfy condition in Eq. (fT8]l . However, for longer loops (as the one shown in Fig. [8]), 
the propagation of these constraints by the Dual algorithm is not perfect, since the information 
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FIG. 7. Energy as a function of the inverse temperature /3 for a 64 x 64 2D EA model, with both types of 
interactions, Gaussian and bimodal. Fuh hnes represent the exact thermodynamical energy as obtained by 
a Monte Carlo simulation, points are the energies obtained under the Dual approximation, and dashed lines 
are the average case energies. 



degrades with distance beyond the plaquette level. Then we need to enforce the constraints of 
Eq. (jlSp by a proper algorithm. At each stage of the freezing process, we define the clusters of 
frozen links as follows: if two frozen links share a spin, then they belong to the same cluster. In 
Figure [8] a cluster of frozen links is represented by bold lines. Notice that, once a spin is fixed in 
a cluster, all other spins are fixed as well by the frozen correlations. On the other hand, different 
clusters of spins can have arbitrary relative orientations. 



FIG. 8. Even if the link marked by the arrow is not the most polarized link according to the marginals 
provided by the Dual algorithm, the spins it connects are fully correlated by the fact that they belong to a 
cluster of frozen links (bold lines). Therefore, the marked link must be immediately fixed accordingly. 

Consider now the situation depicted in Figure[8]and focus on the value of the correlation between 
the two spins connected by the link marked by an arrow. From the fact these two spins belong 
to the same cluster of frozen links (shown as bold link in Figure [8]) we know they are perfectly 
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correlated, however by running the Dual algorithm we could get a weak value for this correlation 
and then proceed by freezing a different link. A set of sub-optimal choices of this kind may finally 
produce a configuration of frozen links where the constraints in Eq. (jlSp are not all satisfied. In 
order to avoid these constraint violations we force any link whose spins are already part of the 
same cluster to be polarized accordingly. The freezing algorithm, therefore, works as follows. 
1: repeat 

2: run the Dual algorithm until convergence (at a low enough temperature) 

3: find the link L with largest finite Jl = Jl + U-p-^ l + Uc-^ l 

4: freeze that link by setting Jl ^ Sign( J^) oo 

5: if link L is connected to clusters C and C of frozen links then 

6: merge clusters C, C and link L in a unique cluster 

7: else 

8: if link L is connected to a single cluster C of frozen links then 

9: add link L to cluster C 

10: else 

11: create a new cluster with link L 

12: end if 

13: end if 

14: for all non-frozen links L' at the boundaries of a cluster of frozen links do 

15: if link L' shares both spins with the same cluster then 

16: freeze link L' accordingly {To avoid violations of constraint} 

17: end if 

18: end for 

19: until all links are frozen 

20: return the spins configuration obtained by setting one spin and fixing the rest according to 
frozen links 

The results obtained with this freezing procedure are quite good. In figure [9] we compare the 
resulting ground state energies with the exact solutions obtained using a web service running an 
exact solving algorithm [l7|. We used an ensemble of 100 EA models on the 2D square lattice 
with Gaussian interactions (so the ground state is not degenerate) and with bimodal interactions. 
Most points are along the bisecting line, meaning that the ground state found by either methods 
are the same. The relative error for the ground state energy is 0.0013 for Gaussian systems, and 
0.00078 for bimodal systems. Looking at how many links are frustrated in one of the solutions and 
unfrustrated in the other, we found that the Dual-|-Freezing algorithm returns 94.3% of the correct 
link correlation signs with respect to the true ground state solution for the Gaussian models. For 
the bimodal system, given the degeneracy of the ground state configuration, it is more probable 
to find the actual ground state energy but for the same reason the link overlap between the exact 
Ground State and the one found with Dual-|- Freezing is significantly lower (86.0%). 

In general we believe these results on the ground states to be very encouraging, considering 
that the Dual algorithm is very fast and not restricted to the 2D case (at variance to fast exact 
algorithms for computing ground states). They provide evidence that the marginals obtained by 
this Dual algorithm are reliable even at very low temperatures. 



V. GENERALIZATION TO OTHER DIMENSIONS 



Let us now consider the region graph based approximation to the free energy for a generic 
D-dimensional (hyper-)cubic lattice, using the same hierarchy of regions: square plaquettes, links 
and spins. After computing the counting numbers for a general D-dimensional lattice, see Eq. ([3]), 
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FIG. 9. Correlation of the Dual+Freezing ground state energy with the exact ground state energy in 
A'' = 16x 16 systems. The top left points correspond to 100 bimodal systems Jjj = ±1, while the right bottom 
points correspond to 100 systems with Gaussian interactions. For bimodal interactions, the degeneracy 
of the ground state improves the probability of actually finding the correct ground state energy (93%), 
and conversely reduces the expected link correlation overlap with the exact ground state solution (86%). 
For Gaussian interactions, the ground state is not degenerated, and only in ~ 7% of the samples the 
Dual+Freezing method finds the actual ground state; however the average link overlap is very high (94%). 
The line f{x) = x is shown to guide the eye. Kindly note that two set of axes are being used. 



the free energy approximation becomes 

-(2L>-3)^5]6i(ai)log 



exp(-^£;p((7p)) 



Plaquettes 



+(2^2 _ 4^ + 1) ^ ^ ^.(g.) log 



exp(-/3EL(aL)) 
hi{si) 



Links 



(19) 



e-x.^{-l3Ei{si)) 



Spins 



Plaquettes are still the biggest regions considered at so have counting number 1, but now each link 
is contained in 2{D — 1) plaquettes, and each spin is in 2D links and 2D{D — 1) plaquettes. The 
message passing equations for the Dual algorithm in D dimensions are then 



U-p-^L = arctanh 



'2(D-1)-1 

tanhyS I Uui^u + Ju 

i 

^2(D-1)-1 \ /2(D-1)-1 

tanh/3 I C/7^,^i^ + Ji^ tanh/3 ^ Uvi^D + Jn 



, (20) 
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where Ui (resp. TZi and Vi) are the 2{D — 1) — 1 plaquettes containing the hnk U (resp. R and D) 
excluding plaquette V. 

In the high temperature phase, this Dual approximation with all u = should be still a valid 
approach for any dimensionality D. At low temperatures, however, the EA model in more than two 
dimensions have a spin glass phase transition and, therefore, we expect the Dual approximation to 
become poorer, as it can not account for a non trivial order parameter. 

By running the Dual algorithm for the 3D EA model we have found a divergence of [/-fields 
around (3 ~ 0.39 for bimodal couplings and around (5 ~ 0.41 for Gaussian couplings. This diver- 
gence is due to the fact the fZ-fields get too much self-reinforced under iteration. This divergence 
does not come as a surprise given that it happens also when studying the simpler pure ferromag- 
netic Ising model. However in the ferromagnetic model the temperature at which [/-fields diverge 
is always below the critical temperature and so the Dual algorithm still provides a very good 
description of the entire paramagnetic phase. 

Unfortunately in the 3D EA the divergence of [/-fields takes place well above the critical tem- 
perature (which is Tc ~ 1.12 for bimodal coupling and at Tc — 0.95 for Gaussian couplings, see 
Ref. for a summary of critical temperatures in 3D spin glasses) and this would make the Dual 
algorithm of very little use. We have studied the origin of this divergence and we have found a 
general principle for reducing the divergence of [/-fields due to self-reinforcement, thus improving 
the convergence properties of the Dual algorithm. The idea is the following. When writing the 
Dual approximation as a constraint satisfaction problem with a non uniform prior, see Eq. ()15p . the 
constraints may be redundant. This is the case for the 3D cubic lattice: indeed both the number 
of links (i.e. variables in the dual problem) and the number of plaquettes (i.e. constraints in the 
dual problem) are 3A^. So, if constraints were independent, the entropy would be null at /? = 
and negative for /3 > (and this is clearly absurd). The solution to the apparent paradox is that 
constraints are not independent: actually only 2/3 of these are independent, and the remaining 
third is uniquely fixed by the value of the former. In this way the correct entropy is recovered 
at /? = 0, given that a problem with 3A^ unbiased binary variables subject to 2A^ independent 
parity-check constraints has entropy A^log(2). The dependence among constraints can be easily 
appreciated by looking at the 6 plaquette around a cube: if 5 of the 6 constraints are satisfied, 
then the sixth one is automatically satisfied and redundant. 

The general rule for improving the convergence of the Dual algorithm is to remove redundant 
constraints (thisprinciple is similar to the maxent-normal property of region based free energy 
approximations [6(]). Redundant constraints have no role in determining the fixed point values for 
the beliefs (since they are redundant), but during the iterations they provide larger fluctuations to 
messages and may be responsible for the lack of convergence. In practice, on a 3D cubic lattice, we 
may remove redundant constraints in many different ways: the basic rule states that one constraint 
(i.e. a plaquette) should be removed for each elementary cube, otherwise if a cube remains with its 6 
plaquettes at least one redundant constraint will exist. We have used two different ways of removing 
one constraint per cube and we have found the same results in the entire paramagnetic phase. So, 
for simplicity, we are going to present data obtained by removing all constraints corresponding to 
plaquettes in the xy plane. 

The Dual algorithm for the 3D EA model on the cubic lattice with no redundant constraints 
converges for any temperature above T ~ 0.8 and so we can use it to study the entire paramagnetic 
phase. The lack of convergence deep in the spin glass phase is to be expected. Just as in the 2D 
case, the Dual algorithm (when converges) still finds the same solution obtained by a Double Loop 
algorithm, and again, it finds the solution nearly 100 times faster (see Fig. [T0|) . Double Loop has 
the apparent advantage of converging at any temperature even at very low ones. However, deep 
in the spin glass phase, where the underlying paramagnetic approximation is clearly inaccurate, 
we believe that an algorithm (like the Dual one) that stops converging, is providing an important 
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FIG. 10. Running times of the Double Loop algorithm 0, [3 (libdai) and of the Dual algorithm on a 
8x8x8 EA model with bimodal interactions [Jij = ±1). The Dual algorithm is generally several orders 
of magnitude faster and returns the same solution as the Double Loop algorithm. 



warning that something wrong is probably happening. Such a warning would be lacking by using 
a Double Loop algorithm. 

In Figure [11] the correlations predicted by the Dual approximation, and those obtained by a 
Parallel Tempering Monte Carlo simulation are compared. At high temperatures the correspon- 
dence is quite good, but not as good as in 2D. However, it is important to stress that the 3D EA 
model is much more difficult to simulate than the 2D case: there is no exact method for com- 
puting the thermodynamics (at variance to the 2D case) and Monte Carlo methods require huge 
thermalization times, while the Dual algorithm runs in linear time with the system size. 

In Figure [12] we show the estimates for the energy obtained from the Monte Carlo and the 
Dual algorithm (both on a single sample and on the average case). The very strong agreement 
between the Dual algorithm results on single samples and on the average case is telling us that 
f7- messages arriving at a given point on the lattice are uncorrelated to a very large extent. In 
other words, the effect of short loops in the lattice is not manifestly present in correlations between 
messages. On the contrary, the comparison between Dual algorithm results and Monte Carlo results 
is good only at high temperatures, and it degrades when approaching the critical temperature. This 
discrepancy can be understood as due to a growing correlation length in the EA model that diverges 
at the critical temperature: our Dual approximation does not account for correlations beyond the 
plaquette level and so it becomes inevitably poorer when the correlation length diverges. However, 
given the extremely fast converging times of the Dual algorithm, it can be viewed as a very effective 
algorithm for sampling the high temperature paramagnetic phase and as a reasonable approximate 
algorithm when approaching the critical point. 
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FIG. 11. Comparison between the correlations (siSj)Duai obtained with the Dual algorithm and the (nearly) 
exact correlations (sjSj )pt obtained with a Parallel Tempering simulation in a 3D EA model of size 8x8x8 
with random bimodal interactions Jij = ±1. At each temperature the correlation coefficient p is reported. 
For the lowest temperature shown, j3 = 1.0, we also report the fraction w-s of pairs of spins such that 

(SiSi)Dual(SiSi)pT < 0. 



VI. CONCLUSIONS 



Wc have introduced a novel Dual algorithm to compute marginals probabilities in the param- 
agnetic phase of frustrated spin models (e.g. spin glasses) on finite dimensional lattices. Inspired 
by the fact that in a paramagnetic phase with no external field each variable is unbiased (i.e. local 
magnetizations are null), the Dual algorithm is derived by adding such paramagnetic constraints 
in the GBP equations. While BP (i.e. Bcthe approximation) and GBP algorithms have serious 
convergence problems at low temperatures even in the paramagnetic phase, the Dual algorithm 
converges very fast in a much wider range thanks to these constraints. The Dual algorithm can 
also be seen as BP on the dual lattice, where the interactions Jij act as external fields on dual 
variables, thus improving convergence properties of the message passing algorithm. 

We have tested the Dual algorithm for the Edwards Anderson spin glass model with bimodal 
and Gaussian couplings on 2D (square) and 3D (cubic) lattices. The results are very encouraging, 
showing convergence in the whole paramagnetic phase (and even slightly in the frozen phase for 
the 3D EA model) and comparing very well with exact correlations measured in Monte Carlo 
simulations. A comparison with a Double Loop algorithm (which is the state-of-the-art among 
general purpose inference algorithms) shows that both algorithms found the same result, but our 
Dual algorithm runs roughly 100 times faster. We also tried to push the Dual approximation to 
the limit, and we used the correlations inferred by the Dual algorithm to compute ground states 
configurations in the 2D EA model by a freezing procedure. Again, we showed that the ground 
states obtained in this way compare very well with exact computations. 

The success of our proposal clearly shows that as long as variables are not long range correlated. 
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FIG. 12. The energy predicted by tlie Dual approximation in 3D EA model, compared to the average case 
energy, and the Monte Carlo simulation. We used a 8 x 8 x 8 system with both types of random interactions, 
bimodal (Jy — ±1) and Gaussian distributed. 

the computation of correlations in a generic spin model can be done in a very fast way by means 
of message passing algorithms, based on mean-field like approximations. This kind of inference 
algorithms do not provide in general an exact answer (unless one uses it at very high temperatures 
or on locally tree-like topologies), and so they can not be seen as substitutes for a Monte Carlo (MC) 
sampling. However there are many situations where a fast and approximate answer is required more 
that a slow and exact answer. Let us just make a couple of examples of these situations. On the 
one side, if one need to sample from very noisy data, an approximated inference algorithm whose 
level of approximation is smaller than data uncertainty is as valid as a perfect MC sampler. On 
the other side, if one need to use the inferred correlations as input for a second algorithm (as for 
the freezing algorithm in Section IIV Bp that will eventually modify/correct these correlations, a 
fast and reasonably good inference is enough. 

The promising results shown in the present work naturally ask for an improvement in several 
directions. For example, in the paramagnetic phase of a model defined on a 3D lattice, our inference 
algorithm could be improved by using the 2x2x2 cube as the elementary region, instead of the 
plaquette. An even more important improvement would be to extend the applicability range of the 
algorithm to the low temperature phase: but this requires a rather non trivial modification, since 
in low temperatures phase the assumption of zero local magnetizations needs to be broken. 
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